/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal.hxx"

#define LINAL_ERROR_TOL 1.0e-9

int main(int argc, char **argv) {


  int uplo, diag;
  double sample, norm, inv_norm_est, inv_norm, diff;
  linal::Flat_ A, B, C, p;
  
  uplo = FLA_LOWER_TRIANGULAR;
  diag = FLA_UNIT_DIAG;

  if (argc != 2) {
    printf("Try :: norm [N]\n");
    return -1;
  }

  int N = atoi( (argv[1]) );
  int datatype = LINAL_DOUBLE_REAL;

  // ---------------------------------------
  // ** Initialization  
  FLA_Init();


  // ---------------------------------------
  // ** Matrices

  // test sample
  A.create    (datatype, N, N);
  p.create    (LINAL_INT,  N, 1);

  // inverse matrix
  B.create    (datatype, N, N);
  C.create    (datatype, N, N);

  // random setup and factorization
  FLA_Random_matrix(~A);
  FLA_LU_piv(~A,~p);
  FLA_LU_nopiv(~A);

  // get inverse
  FLA_Copy(~A, ~B);
  FLA_Triangularize(uplo, diag, ~B);
  FLA_Set(FLA_ZERO, ~C);
  for (int i=0;i<N;++i)
    C(i,i)= 1.0;

  FLA_Trsm( FLA_LEFT, uplo, FLA_NO_TRANSPOSE,
	    diag, FLA_ONE, ~B, ~C );

  // ---------------------------------------
  // ** Check

  // calculate norm
  norm         = linal::norm1(uplo, diag, A);
  inv_norm_est = linal::inv_norm1(uplo, diag, A);
  inv_norm     = linal::norm1(uplo, diag, C);
  diff         = inv_norm_est - inv_norm;
  
  int rval;
  printf("- TEST::");
  for (int i=0;i<argc;++i)
    printf(" %s ", argv[i] );
  printf("\n");

  if (diff < LINAL_ERROR_TOL) {
    printf("PASS::Norm %E, Inv_Norm_Est %E, Inv_norm %E, Diff %E ( %lf percents error)\n",
	   norm, inv_norm_est, inv_norm, diff, -1.0*diff/inv_norm*100);
    rval = 0;
  } else {
    printf("FAIL::Norm %E, Inv_Norm_Est %E, Inv_norm %E, Diff %E ( %lf percents error)\n",
	   norm, inv_norm_est, inv_norm, diff, -1.0*diff/inv_norm*100);
    rval = -1;
  }

  // ---------------------------------------
  // ** Matrices
  C.free();
  B.free();

  p.free();
  A.free();

  // ---------------------------------------
  // ** Finalization
  FLA_Finalize();

  return -1;
}
